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Abstract 


The present paper will present a description of and results from a new solution 
algorithm for the compressible Navier-Stokes equations developed by a team of re- 
searchers at Langley Research Center. The main features of the algorithm are second- 
or third-order accurate upwind discretization of the convective and pressure deriva- 
tives and a relaxation scheme for the unfactored implicit backward Euler time method, 
implemented in a finite-volume formulation. 

Upwind methods have been successfully used to obtain solutions to the Euler 
equations for flows with strong shock waves. One reason for this success is that 
these methods directly simulate the signal propagation features of hyperbolic 
equations. Furthermore, these methods have the advantage of being naturally dis- 
sipative. Upwind differencing has been used for the pressure and convective terms in 
the present Navier-Stokes algorithm while central differencing has been used for the 
viscous terms. The particular upwind method being used is based on the flux-vector- 
splitting technique developed by Van Leer (ref. 1) and both second- and third-order 
accurate discretizations have been developed. 

Currently, the most widely used implicit solution techniques for the Navier- 
Stokes equations use approximate factorization (AF) methods to treat multi- 
dimensional problems. Although the implicit AF methods are unconditionally stable in 
two dimensions, they require an optimal set of iteration parameters, which are diffi- 
cult to obtain, for rapid convergence. Furthermore, the implicit AF schemes are only 
conditionally stable for the three-dimensional Navier-Stokes equations. However, the 
upwind discretization leads to a diagonally dominant matrix structure which allows 
large time steps to be taken in multi-dimensional problems. The time integration 
scheme being used in the present algorithm corresponds to a line Gauss-Seidel re- 
laxation method. This method produces good convergence rates for steady-state flows, 
and most of the algorithm has been vectorized on the NASA Langley VPS 32 computer. 

The Navier-Stokes algorithm has been tested for several two-dimensional flow 
problems such as the laminar boundary layer flow on a flat plate and the flow pro- 
duced by an oblique shock wave impinging on a laminar boundary layer developing on a 
flat plate. Solutions for both problems gave excellent results which are presented 
in this paper. Present effort is directed toward the extension of the scheme to the 
full three-dimensional Navier-Stokes equations. 
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Upwind Relaxation Algorithms for the 
Euler/Navier-Stokes Equations 

The purpose of this paper is to present upwind relaxation algorithms for accur- 
ate and reliable steady-state solutions to either the compressible Euler (inviscid) 
or Navier-Stokes (viscous) equations. The main feature of the algorithms is the use 
of second- or third-order accurate upwind discretizations of the convective and pres- 
sure terms, which enables implicit line relaxation strategies to be developed for 
rapid solutions to the steady-state equations. The basic algorithm for three- 
dimensional flows is discussed. Both inviscid and viscous (using thin-layer Navier- 
Stokes approximations) computations are shown for a series of two-dimensional flows. 


91 



3-D Relaxation Algorithm 


The relaxation algorithm for the time-dependent Euler equations written in con- 
servation form and generalized coordinates is shown. The equation to be solved with 
relaxation corresponds to the linearized, backward-time approximation in delta form. 
The inviscid flux is split according to the eigenvalues of the characteristic equa- 
tion, which recognizes the signal propagation features of the hyperbolic equations 
and enables relaxation approaches to be exploited for solutions to the steady-state 
equations. The particular flux splitting technique used corresponds to that devel- 
oped by Van Leer (ref. 1), although the advantages resulting from upwind differencing 
would apply equally as well to other splitting techniques. Applying upwind differ- 
encing in the streamwise ( £) direction on the implicit (left-hand) side of the 
equation in delta form leaves an equation to be solved in the cross flow (tvC) 
plane. The streamwise relaxation is effected by sweeping in the ^-direction 
through the mesh, alternating the direction of the sweep every other pass in order to 
maintain stability for higher order differencing. 

Non-linear updating of the residual is indicated, corresponding to upstream 
cross-flow planes evaluated at time level n+1 when sweeping from upstream to down- 
stream, and vice-versa. The upwind relaxation scheme indicated is unconditionally 
stable and maximal damping occurs at large time steps. This is in contrast to 
approximately factored approaches, which require an optimal set of iteration 
parameters, which are difficult and time consuming to obtain, for rapid 
convergence. The relaxation scheme also has the advantage that it recovers 
conventional space marching techniques for supersonic flows in the ^-direction. 
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Solution of Cross-Flow Equations 


The upwind differenced equation in the cross-flow plane can be solved using 
either relaxation or approximate factorization. The upwind relaxation shown corre- 
sponds to alternate line Gauss-Seidel sweeping with nonlinear updating of the re- 
sidual indicated, although in practice the residual is updated linearly using the 
Jacobian matrices from the implicit side of the equation. With first-order implicit 
upwind differencing, either approach leads to two sweeps across the plane, solving a 
system of tridiagonal equations on each line. With relaxation in the cross-flow 
plane, the scheme is unconditionally stable for the linear wave equation and maximum 
damping occurs with large time steps. However, because the algorithm is recursive in 
nature, it cannot be completely vectorized on current pipeline supercomputers such as 
the CDC CYBER 205. The time step is limited to obtain optimal convergence with the 
approximate factorization in the cross-flow plane, but the scheme is completely 
vectorizable. Thus, a tradeoff exists between convergence rate and computational 
rate in determining the most efficient strategy. Vectorizable relaxation strategies, 
such as checkerboard algorithms, lead to faster computational rates with slightly 
lower convergence rates, owing to the underrelaxation required for stability. 
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Spatial Differencing 


The convective and pressure terms are split according to the flux vector split- 
ting technique developed by Van Leer, which is implemented in a control- volume form- 
ulation, The flux difference is computed from the difference of split fluxes across 
all boundaries, where the upwinding is incorporated through an interpolation of 
conserved variables. A one-parameter family of interpolations is used. This family 
ranges from fully upwind second-order to upwind-biased third-order approximations. 
Viscous terms are differenced using second-order central differences, so that the 
method for viscous flows is limited currently to second order. 


Upwind flux splitting for convective and pressure terms 
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Inviscid Shock Reflection From Flat Plate 


Application of the relaxation algorithm to the solution of an inviscid shock 
reflection from a flat surface is shown. The freestream Mach number is 2.9 and the 
shock angle is 29°. Fully upwind spatial differencing in the streamwise direction is 
used along with vertical Gauss-Seidel line relaxation? this leads to the quadratic 
convergence shown since the flow is fully supersonic. First-order differencing leads 
to block tridiagonal line inversions; second- or third-order differencing leads to 
block pentadiagonal line inversions. The residual history shown is for a global line 
relaxation strategy (sweeping downstream through the mesh). Further improvements in 
efficiency can be realized by removing linearization errors completely on a line 
before moving to the next line, and, thus, recovering a space marching method. The 
space marching strategy has also been extended to efficiently treat supersonic flows 
with embedded subsonic regions. 
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Pressure Contours For Inviscid Shock Reflection 


Pressure contours for the inviscid shock reflection are shown with various 
orders of spatial differencing. The grid for each case is uniformly spaced in each 
of the two directions. The higher order methods are much less dissipative then the 
first-order scheme and show much higher resolution of the shock wave. 
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Residual Histories for First- and Second-Order Schemes 


Residual histories for the transonic flow over a bump in a channel are shown. 

The incoming Mach number is 0.85, and a significant region of transonic flow exists 
as shown subsequently. The particular results shown were obtained with approximate 
factorization (AF) and with vertical line relaxation (LR) methods on a 34 x 18 H-mesh 
with uniform spacing. The Courant number for the LR scheme is on the order of 200, 
while that of the AF scheme is limited to 20. One iteration refers to two passes 
through the mesh using AF and one pass through the mesh using LR. The computational 
work per iteration on a scalar processor is roughly equal for the two methods and 
thus the LR scheme is much more efficient. On a vector processor, the computational 
rate per iteration for the LR scheme is two to three times slower than the AF scheme, 
even though more than 90 percent of the LR scheme has been vectorized. The only 
scalar operations in the LR scheme are the recursive forward-backward substitutions 
on a line; the lower-upper decomposition of the matrix equations is completely 
vectorizable. Thus, the relative efficiency of the two schemes depends on the type 
of computer on which the two approaches are implemented. 
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Computational Grid for Bump in Channel 


A computational grid for the transonic flow over a bump in a channel is shown. 
The thickness of the bump is 4.2 percent of the chord, and the channel is approxi- 
mately two chords in height. The computational grid extends two chords ahead of and 
behind the bump and is clustered in the region near the bump. The grid is uniform in 
the vertical direction. 


h/c = 2.073; t/c = 4.2 %; 85 x 41 H-mesh 
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Pressure Contours for Bump in Channel 

Pressure contours for the transonic flow over a bump in a channel obtained on 
the 85 x 41 computational grid are shown. A substantial region of transonic flow 
exists over the bump, and the shock is captured with two transition zones using the 
present method. The results shown were obtained with the fully upwind second-order 
scheme ( k = -1 ) . 


Mqq = 0.85; t/c = 4.2%; h/c = 2.073; 85x41H-mesh 
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One-Dimensional Couette Flow 

Applications of the relaxation algorithm to a series of viscous model problems 
are shown. The first is one-dimensional Couette flow between a fixed lower plate and 
a moving upper plate at low Reynolds number and Mach number. Velocity profiles ob- 
tained with third-order upwind-biased differencing for the convective and pressure 
terms and second-order central differencing for the viscous shear terms are shown 
where the parameter B represents an applied pressure gradient term. No-slip velocity 
and fixed-wall temperature boundary conditions were used. The relaxation method at 
large time steps in one-dimension corresponds to Newton’s method for finding 
solutions to nonlinear systems of equations and convergence of the residual to 
machine zero can be obtained in less than ten iterations. The computed velocity 
profiles agree very closely with the exact result. 
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Flat-Plate Boundary Layer: Second-Order Differencing 

The accurate computation of viscous effects with an upwind method is an import- 
ant concern. Results are shown for the laminar flow over a flat plate using second- 
order (k= -1) differencing for the pressure and convective terms. The freestream 
Mach number is 0.5 and the Reynolds number based on the length of the plate is 
10,000. The results were obtained on a 30 x 40 computational grid, which is uniform 
in the streamwise direction and stretched in the direction normal to the plate. The 
computational domain is sketched; inflow conditions were applied upstream of the 
plate and outflow conditions above and downstream of the plate. No-slip, adiabatic 
wall conditions were specified on the plate. The computations recover the similarity 
velocity profiles, which develop very quickly downstream of the leading edge. The 
second-order scheme resolves the boundary layer adequately with approximately twenty 
points in the layer. 
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Flat “Plate Boundary Layer: Third-Order Differencing 

Results are shown using third-order upwind-biased differencing (k = 1/3) for the 
convective and pressure terms. The accuracy relative to the second-order differenc- 
ing is evident by comparison with the previous figure. The third-order differencing 
can be implemented with only a very small additional computational effort relative to 
the second-order scheme, since the computational molecule for both schemes only in- 
volves five points in any one direction. Central difference methods and the present 
upwind method are formally of the same order (second order) for viscous flows and one 
should expect roughly the same truncation error level. This has been verified in the 
present effort both analytically and numerically. In contrast, first-order upwind 
differencing leads to a boundary layer thickness four to five times that shown. 
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Shock-Boundary Layer Interaction 


The final problem is that of an oblique shock wave impinging on a laminar bound- 
ary layer developing on a flat plate. The conditions correspond to the experiments 
of Hakkinen et al. (ref. 2) at a free-stream Mach number of 2.00 and a Reynolds 
number based on the length from the leading edge to the shock impingement point of 
2.96 x 10 . The shock is of sufficient strength to cause a separation of the de- 
veloping laminar boundary layer, as evident in the computational results shown using 
thin-layer approximations to the Navier-Stokes equations. The pressure contours show 
the compression and expansion waves downstream of the shock impingement, which are 
caused by the interaction of the incoming shock wave with the separated boundary 
layer. 
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Shock-Boundary Layer Interaction: Grid Convergence 


In order to assess the accuracy of the present upwind method, results from a 
grid convergence study are shown. Mesh densities two and three times the original 
grid were used. Details of the separated flow region are generally resolved with the 
second mesh, although a small difference exists between the second and third meshes 
near reattachment. Outside the separated flow region, very little difference in wall 
pressure exists in the results shown. The details in the present calculations 
clearly show the pressure plateau in the separation region. 
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Shock-Boundary Layer Interaction: Skin Friction Coefficient 


The skin friction coefficient predicted with the present method for the three 
meshes is shown. The skin friction exhibits three local extremum, two minimums and 
one maximum, in the separation zone. The separation extends from approximately 0.75 
to 1.20 reference lengths from the leading edge of the plate. 
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Shock-Boundary Layer Interaction: Wall Pressure and Skin Friction 

A comparison of wall pressure and skin friction from the present results with 
experiment is shown. The pressure plateau and the extent of separation are closely 
predicted. The arrows in the figure indicate that only the presence of separation 
was deduced from the experiment; the actual level of skin friction in the separation 
was not measured. 
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Concluding Remarks 


Efficient solutions to steady-state Euler and Navier-Stokes equations are pos- 
sible through the upwind relaxation algorithms discussed previously. The methods 
recover conventional space marching methods for fully supersonic flows and can 
be shown to be unconditionally stable in three dimensions. Accurate solutions to 
several inviscid and viscous model problems, including flat -plate boundary layer flow 
and shock/boundary layer interaction with separation, have been demonstrated. The 
method is currently being implemented in three dimensions. 


Efficient solutions to steady-state equations possible 
through upwind relaxation 

Recovers space marching method for supersonic flows 
Unconditionally stable in three dimensions 

Accurate solutions to several viscous flows demonstrated 
using second/third-order upwind differencing for 
convective and pressure terms 

Flat plate boundary layer 

Shock/boundary-layer interaction with separation 
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